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We present a new phenomenological gravitational waveform model for the inspiral and coales- 
cence of non-precessing spinning black hole binaries. Our approach is based on a frequency domain 
matching of post-Newtonian inspiral waveforms with numerical relativity based binary black hole 
coalescence waveforms. We quantify the various possible sources of systematic errors that arise 
in matching post-Newtonian and numerical relativity waveforms, and we use a matching criteria 
based on minimizing these errors; we find that the dominant source of errors are those in the 
post-Newtonian waveforms near the merger. An analytical formula for the dominant mode of the 
gravitational radiation of non-precessing black hole binaries is presented that captures the phe- 
nomenology of the hybrid waveforms. Its implementation in the current searches for gravitational 
waves should allow cross-checks of other inspiral-merger-ringdown waveform families and improve 
the reach of gravitational wave searches. 
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I. INTRODUCTION 

As a generalization of the classic Kepler problem in 
Newtonian gravity, the binary black hole (BBH) system 
in general relativity is of great interest from a fundamen- 
tal physics viewpoint. Equally importantly, this system 
has received a great deal of attention for its relevance 
in astrophysics and, in particular, as one of the most 
promising sources of detectable gravitational radiation 
for the present and future generations of gravitational- 
wave detectors, such as LIGO [1], Virgo '2\, GEO600 "S], 
LISA i4j or the Einstein Telescope ^. The Kepler prob- 
lem can be solved exactly in Newtonian gravity and it 
leads to the well-known elliptical orbits when the sys- 
tem is gravitationally bound. In contrast, in general 
relativity, closed orbits do not exist and the BBH sys- 
tem emits gravitational waves (GWs) which carry away 
energy, thereby causing the black holes to inspiral in- 
wards, and to eventually coalesce. The emitted GWs 
are expected to carry important information about this 
process, and it is one of the goals of gravitational wave 
astronomy to detect these signals and decode them. 

No analytic solutions of Einstein's equations of gen- 
eral relativity are known for the full inspiral and merger 
of two black holes. Post-Newtonian (PN) methods can 
be used to calculate an accurate approximation to the 
early inspiral phase, using an expansion in powers of v/c 
(where v is the orbital velocity and c is the speed of 
light). As for the coalescence phase, starting with |S{8], 
the late inspiral and merger has been calculated by large- 
scale numerical solutions of the full Einstein field equa- 



tions. Since the initial breakthroughs in 2005, there has 
been dramatic progress in numerical relativity (NR) sim- 
ulations for GW astronomy, including many more orbits 
before merger, greater accuracy and a growing sampling 
of the black hole-binary parameter space. A summary of 
the published "long" waveforms is given in the review [S] , 
and a complete catalog of waveforms is being compiled 
at |10j : more recent work is summarized in [llj . NR 
results are now accurate enough for GW astronomy ap- 
plications over the next few years |12j . and have started 
playing a role in GW searches [13l [T4j . 



Given PN and NR results, it is promising to try and 
combine them to produce "complete" inspiral-merger- 
ringdown waveforms. PN techniques in their standard 
formulation become less accurate as the binary shrinks, 
and the approximation breaks down completely some- 
where prior to the merger. NR waveforms, on the other 
hand, become more and more computationally expensive 
the larger the number of cycles that one wishes to sim- 
ulate; the longest published data spans 16 orbits for the 
equal-mass non-spinning case |15j . We therefore would 
hope to combine PN and NR results in the region between 
the point where NR simulations start and where PN 
breaks down. To do this it is critical to verify that the PN 
and NR results are in good agreement in this region and 
that there is a consistent PN-NR matching procedure. 
Much work has been done in comparing PN and NR re- 
sults over the last 5-15 orbits before merger for a variety 
of physical configurations, such as the equal-mass non- 
spinning case p^H22] , the equal-mass non-precessing-spin 
case [23], and the unequal-mass spinning case [21] ■ The 



consistency of PN amplitudes during the merger and 
ringdown has also been studied [HI [25l [26] . These stud- 
ies suggest that a sufficiently accurate combination of PN 
and NR results should be possible. One topic that has 
not received much attention, however, is the systematic 
errors that arc introduced by different choices of match- 
ing procedure. 

One of the aims of this paper is to further understand 
and quantify the various systematic errors that arise in 
the matching procedure. There are thus far two kinds 
of approaches to the PN-NR matching problem, both of 
which have yielded successful results. The first is the 
Effective-One-Body (EOB) approach p7H5Dj. Originally 
motivated by similar techniques in quantum field theory, 
the idea is to map the two body problem into an effec- 
tive one-body system with an appropriate potential and 
with the same energy levels as the two-body system. It 
was shown |27j that the appropriate one-body problem 
(for non-spinning black holes) is that of a single particle 
moving in a deformed Schwarzschild spacetime. It turns 
out that most parameters of this one-body system can 
be found by using the appropriate PN calculations, and 
the remaining parameters are calculated by calibrating 
to NR simulations. This approach has been successful so 
far for non-spinning systems where only a few parameters 
need to be calibrated by NR [31-35 . The spinning case 
is more complicated, and work is underway to extend the 
parameter space described by the model [36j. 

A complementary approach is to perform a phe- 
nomenological matching of the GW waveforms in a win- 
dow (which could be either in the time or frequency do- 
main) where both PN and NR are expected to be good 
approximations to the true waveform. The first step is 
to construct a hybrid PN-NR waveform by matching the 
two waveforms within the matching window. The wave- 
form is completely PN before this window, completely 
NR afterward, and it interpolates between the two in 
the matching window. Once the hybrid waveform is con- 
structed and we are confident about the matching proce- 
dure, the second step is to fit the hybrid waveform to a 
parametrized model containing a number of phenomeno- 
logical coefficients and finally to map them to the physi- 
cal parameters of the system. The resulting model would 
thus be parametrized by the masses and spins of the two 
black holes (and eccentricity if appropriate) . Most of the 
work in this approach has thus far been based on match- 
ing PN and NR waveforms in the time domain, but then 
producing a phenomenological model in the frequency 
domain, which is often more convenient for data-analysis 
applications [57H1D]. See also ^41j for a complementary 
construction. In this paper we take a slightly different 
approach: both the construction of the PN-NR hybrid 
waveform and the matching to a phenomenological model 
are carried out in the frequency domain. The reasons for 
this are twofold. First, we find it easier to work in the 
frequency domain since the quantities used to estimate 
the errors of our matching procedure and the goodness of 
the fit, such as waveform overlaps, are conveniently for- 



mulated in Fourier space. Second, and more importantly, 
in light of the potential errors in the hybrid construction, 
comparing results between two independent methods is 
a valuable way of ensuring that the matching procedure 
is robust. The frequency domain construction presented 
here is complementary to the time domain method of [40j . 

The phenomenological waveform family presented 
in [40, used a simple piecewise ansatz for the phase of 
the hybrid PN-NR waveform, and another for the ampli- 
tude. The resulting analytic model was found to agree 
with the hybrid waveforms with overlaps above 97% for 
most black hole binary systems that would be observ- 
able by the current LIGO detectors. In this paper we 
investigate whether the fidelity of the phenomenologi- 
cal waveforms can be improved by using ansatze that 
make smooth transitions between their inspiral, merger 
and rindown forms. This procedure also allows us to fur- 
ther test the robustness of the phenomenological model's 
construction to variations in its analytic form. 

The main results of this paper are the following. We 
construct hybrid waveforms for binary black hole sys- 
tems with aligned spins in the frequency domain. We 
do this by combining 3.5PN waveforms in the stationary 
phase approximation with a number of NR results. We 
show that this construction is internally consistent and 
it yields hybrids which are, for the most part, sufficiently 
accurate for the initial and advanced LIGO detectors. 
Notably, the difference between the different PN approx- 
imants is a more significant source of error than the nu- 
merical errors in the NR waveforms. Using these hybrid 
waveforms, we construct a phenomenological frequency- 
domain waveform model depending on three parameters 
(as in jini) and covering the space of aligned spins and 
moderate mass ratios. We show that the model fits the 
original hybrid waveforms with the overlaps (maximized 
over the model parameters) better than 97% for Ad- 
vanced LIGO (and for the most part, better than 99%) 
for essentially all black hole systems observable with Ad- 
vanced LIGO, i.e. for systems with total mass ranging up 
to ~ 4OOM0. These results are comparable to those ob- 
tained in [401 , suggesting that the phenomenological con- 
struction is indeed robust and that the phenomenological 
model waveforms are useful for detection purposes. 



Sections [II] and |III| describe the post-Newtonian wave- 
form model and the numerical waveforms that we employ. 
Section [TV] describes the fitting procedure and the vari- 
ous systematic errors that appear in this procedure. It 
quantifies the reliability of the waveforms for specific GW 
detector and signal-to-noise ratios (SNRs). Sectionlvlfits 
these hybrid waveforms to an analytic model. It shows 
that the model provides a good representation of the hy- 
brid waveforms and can be used in GW searches in the 
appropriate parameter space. Finally, section |VI| con- 
cludes with a summary and suggestions for future work. 



II. NUMERICAL SIMULATIONS OF 
NON-PRECESSING BLACK HOLE BINARIES 



found in 



In this section we summarize the numerical waveforms 
used in this paper. Since the first successful numeri- 
cal simulations of equal-mass, non-spinning binary black 
hole mergers were published 0-^ the NR community 
has continued exploring the parameter space of the BBH 
system. Each black hole is described by a mass and a 
spin vector, and the binary's trajectory is described by 
adiabatically evolving Keplerian orbits, so 17 parameters 
are needed to describe the binary system (see e.g. |42|). 
Besides the two masses and spin vectors, we need a fidu- 
cial time t(j and orbital phase 0o at to, the distance to 
the source and its sky-location, two parameters for the 
unit vector normal to the orbital plane, and finally, if 
non-circular orbits are considered, we additionally need 
the eccentricity and the direction of the semi-major axis. 
Sufficiently close to or during the merger, this description 
in terms of Keplerian orbits will break down, and higher 
order black hole multipoles might play a role as well. 

Due to the complexity of this parameter space, most 
match-filtered searches for coalescing binaries have so 
far employed non-spinning templates, neglecting the ef- 
fect of the spin by assuming a small, tolerable loss in 
SNR [48H50] . Dedicated searches for spinning binaries 
have attempted to model an enlarged parameter space 
by using a template family designed to capture the spin- 
induced modulations of the gravitational waveform [5T] . 
In [51] the spin effects were modeled using unphysical 
phenomenological parameters; however, it would be de- 
sirable to devise searches for spinning systems based on 
strictly physical parameters. Indeed, [H] showed that 
from the point of view of detection efficiency at a given 
false-alarm rate, a search based on non-physical spinning 
templates is not superior to a non-spinning search un- 
less specific signal-based vetoes and other tools are de- 
vised. The performance of spinning searches would in- 
crease with the use of templates determined by physi- 
cal rather than phenomenological parameters. That was 
the motivation for the waveform family presented in |40j , 
where as a first step in modeling the full spinning-binary 
parameter space, only binaries with non-precessing spins 
were considered. Additionally, it is known from PN treat- 
ments of the inspiral |53j and from numerical simulations 
of the merger j54] (which though only considers equal 
mass systems) , that the dominant spin effect on the wave- 
form is from the total spin of the system. Indeed, in [40] 
it was found that the effect of the black hole spins can 
be modeled with suficient accuracy using only one spin 
parameter, roughly corresponding to the total spin of the 
two black holes. We adopt the same approach here. 

There are a number of NR simulations of non- 
precessing systems for a variety of spin values and mass 
ratios. Results with the BAM code are reported in [22 for 
the orbital hang-up case and in |55j for anti-aligned spins. 
The CCATIE simulations are presented in [101 liSl F551 157] : 
a long spectral simulation with anti-aligned spins can be 



1. NR waveforms and codes 

The NR waveforms employed in the construction of 
the hybrid model used in this paper are summarized in 
Table |TJ They have been produced with four indepen- 
dent NR codes, BAM, CCATIE, Llama and SpEC. The first 
3 codes use the moving-puncture approach [3 [59] to solve 
the Einstein equations in a decomposed 3-1-1 spacetime 
while the last implements the generalized harmonic for- 
mulation [201 W\- BAM and CCATIE use computational 
domains based on Cartesian coordinates, while the SpEC 
code uses a sophisticated series of spherical and cylindri- 
cal domains; in the wave zone, the outer computational 
domains have the same angular resolution, thus the com- 
putational cost only increases linearly with the radius of 
the outermost shell. A summary of the properties of the 
three codes is given in [T^]. The Llama code jl^l [HT] is 
based on finite differencing but the set-up of the numeri- 
cal grid in the outer wave zone is, as in SpEC, also based 
on spherical coordinates with constant angular separa- 
tion. The large wave-zone enables accurate waveform ex- 
traction at large distances, accurate extraction of higher 
angular modes of the radiation, and it allows the outer 
boundary to be far enough away so that it is causally 
disconnected from the sphere where the radiation is ex- 
tracted. 

The BAM data-set #1 covers the parameter space of 
non-spinning systems for several mass ratios during at 
least the last 5 orbits before merger (length ^ 1100 — 
1450 M, where M is the total ADM mass of the space- 
time) [18l[38l[39l[55]. Data-set #2 consists of moderately 
long simulations covering at least the last 8 orbits before 
merger (length ^ 1500 — 2200 M) for equal-mass systems 
with equal spins, and are described in depth in [^[55] , 
Data-set #3 consists of unequal-mass, unequal-spins sim- 
ulations [40j. Data-set #4 is a simulation with unequal 
mass and unequal spins employed in the verification of 
our fitting mode 0D]. For the sets #1-4, initial momenta 
for quasi-circular orbits were computed for non-spinning 
cases according to the procedures described in [ii] , lead- 
ing to low-eccentricity (e < 0.006) inspiral evolutions. A 
number of different methods were used for the spinning 
cases [23l [55l [62] , depending on which method gave the 
lowest eccentricity for a given configuration. The GW 
radiation is calculated from the Weyl tensor component 
^4 (see e.g. [53]) and extracted at a sphere with radius 
R = 90 M. In all cases the uncertainty in the phase is less 
than 0.1 rad during inspiral (up to Mlo = 0.1), and less 
than 0.5 rad during merger and ringdown. The uncer- 
tainty in the amplitude is less than 0.5% during inspiral, 
and less than 5% during merger and ringdown. 

The CCATIE data-sets #5, #6 and #7ab correspond to 
the S-, U-, r- and i-sequences studied in [54]. They 
span the last ~ 4 — 5 orbits before merger (length 
^ 500 — 1000 Af ) and are in fact not sufficiently long 



TABLE I. NR codes and configurations used for the construction and verification of our fiybrid waveforms and phenomenological 
model. The mass ratio q is defined as mi/m2, assuming m\ ^ 7712; Xi,2 are the dimensionless spins defined in Eq. (3.1 1; a 
positive value of xi,2 means that the spin is aligned with the orbital angular momentum L, and negative values are anti-aligned. 



Data Set 


Code 


Mass ratios 


#1 


BAM gam] 


qe {1,1.5,2,2.5,3 


#2 


)? 


q = l 


#3 


!7 


gG{2,3,4} 


#4 


•>•> 


q = 3 


#5 


CCATIE pS] 


g = l 


#6 


11 


g = l 


#7ab 


11 


q = l 


#8 


Llama [iB] 


ge{i,2} 


#9 


SpEC [47] 


g = l 



Spins 



Extraction of GW signal 



(xi,X2) = (0,0) 

(Xi, X2) = ia,a),ae ±{0.25, 0.5, 0.75, 0.85} 

(xi,X2) = (a,a), a G {±0.5,0.75} 

(Xi,X2) G {(-0.75, 0.75), (0,0.8333)} 

(Xi,X2) = (a,a), a €{0,0.2,0.4,0.6} 

(Xi,X2) = (a,-a), a €{0,0.2,0.4,0.6} 

(Xi, X2) = (±0.6, a),ae {±0.3, 0, -0.6} 

(xi,X2) = (0,0) 

(Xi,X2) = (0,0) 



at r = 90M 



at r = 160M 

■>1 
•>1 

Null InfinitjF] 
at r ^ ocj^ 



^ Only the GW radiation corresponding to the Llamia q = 1 simulation has been extracted at future null-infinity using the 

Cauchy-characteristic method; the q = 2 waveform has been extracted at finite radius and extrapolated to r — > 00. 
^ Using the extrapolation method described in |15| with extrapolation order n = 3. 



for use in the hybrid construction. They are still useful 
to independently verify the reliability of our phenomeno- 
logical fit. Data-set #5 corresponds to the hang- up con- 
figuration analogous to the BAM set #1; data-set #6 con- 
sists of configurations virith (xi,X2) — (fli—o), i-e. zero 
net spin; data-set #7a was analyzed in [45j in the con- 
text of the study of the recoil velocity ("kick") of the 
final merged black hole. GW radiation is extracted at 
R = 160 M via the Regge-Wheeler-ZeriUi formalism for 
perturbations of a Schwarzschild black hole 64-67; . 

Data-set #8 consists of two waveforms for non- 
spinning black holes with mass ratios q — 1,2. The black 
holes are evolved with the Llama code according to the 
set-up reported in [IB]. The outer boundary is placed at 
3600 Af and the initial separation is IIM, corresponding 
to 8 orbits in the inspiral phase followed by merger and 
ringdown. Wave extraction for the q = 1 configuration is 
done via the Cauchy-characteristic method [Ml 112]; tak- 
ing boundary data from the numerical spacetime for a 
subsequent characteristic evolution of the metric to null- 
infinity, thereby obtaining waveforms that are mathemat- 
ically unambiguous and free of any systematic finite ra- 
dius and gauge effects. The only remaining source of 
error is due to numerical discretization. The equal-mass 
waveform produced with this code was reported in [55] . 
while the q = 2 waveform is new. For these data, the 
uncertainty in the phase is comparable to that of the BAM 
waveforms, while the uncertainty in the amplitude is at 
least an order of magnitude lower, because of the more 
sophisticated wave extraction procedure [JBl [551 [SH] ■ 

Data-set #9 consists of a long non-spinning, equal- 
mass simulation that follows 16 orbits of the binary 
plus merger and ringdown of the final black hole (length 
~ 4300 Af). These are publicly available data 70 which 
were originally computed using the SpEC code with neg- 
ligible initial orbital eccentricity ('-^ 5 x 10~^). The GW 



radiation is extracted via ^4 in a similar manner to #1- 
4 and extrapolated to infinity. The phase uncertainty is 
less that 0.006 rad during inspiral, and less than 0.02 rad 
during merger and ringdown; the amplitude uncertainty 
is less than 0.1% during inspiral, and less than 0.3% dur- 
ing merger and ringdown. A full description of this simu- 
lation is given in [15) . The long duration of the waveform 
allows for its use in the estimation of the errors associated 
with the length of the NR data. In particular, since it 
contains physical information at lower frequencies, it can 
be matched to PN results at lower frequencies, where the 
PN errors are expected to be smaller (see the discussion 
around the right panel of Fig. [6]) . 



2. Going from ^4 to h 

The gravitational waveforms calculated using NR 
codes are typically reported in terms of the Weyl ten- 
sor component "^4, which is a complex function that en- 
codes the two polarizations of the outgoing transverse 
radiation. ^4 is related to the two polarizations of the 
gravitational wave perturbation h^^ (in the transverse- 
traceless gauge) via two time derivatives 



*4 



dt' 



[h+{t) - ih^it)] 



(2.1) 



Going from \l/4 to ft.+,x thus involves two time integra- 
tions and requires us to fix two integration constants ap- 
propriately, corresponding to the freedom to add a linear 
function to the strain. 

The frequency domain offers a straightforward way of 
calculating the strain 



ihy 



(2.2) 
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FIG. 1. These figures demonstrate the strain waveform ob- 
tained by the frequency domain division method of calculating 
h from ^4. We consider the NR simulation from data-set #8 
of Table|l]with q — 1, and we start with the dominant mode of 
5'4 from this simulation. The upper panel shows A = r\h{f)\ 
(where r is the extraction radius) obtained by the frequency 
domain division, and the lower panel shows h+(t). The win- 
dow function employed by our inverse Fourier transform al- 
gorithm is responsible for the partial loss of the first cycle of 
the waveform. Nevertheless, a clean \h{t)\ during the rest of 
the inspiral is observed. 



from ^4, since integration is replaced by division: 



lNRfr\ _ ^4 (/) 



^^NR(^)g.*-(/) 



(2.3) 



where x(f) denotes the Fourier transform of x(t) as de- 
fined in Eq. (4.2). In the limit of large signal durations, 



the integration constants only affect the zero-frequency 
component of the signal in the frequency domain. For 
finite duration signals, the effect of the integration con- 
stants will spill over into higher frequencies like 1//. 
Since for our purposes, the numerical simulation provides 
useful information only starting at a finite frequency, we 
conveniently apply a high-pass filter to the data, thus 
reducing the effect of the integration constants without 
using a fitting procedure. When performing the divi- 
sion in the frequency domain, tanh- window functions are 
employed to pass-filter the data before computing the 
Fourier transform. Fig. [T] illustrates the efficacy of our 
approach for the Llama equal-mass waveform. Though 
we do not discuss it further here, in general we find that 
the time- and frequency-domain integration techniques, 
both with some fine-tuning, yield comparable results. 



III. ANALYTICAL WAVEFORMS FOR 

SPINNING BINARIES USING THE 

POST-NEWTONIAN APPROACH 

Coalescing compact binaries such as BBHs can be ac- 
curately modeled by the FN approximation to general 
relativity at least during the major part of the long inspi- 
ral phase, under the assumptions of a weak gravitational 
field [7T]. In order to obtain an analytical description of 
the early inspiral in the Fourier domain we construct the 
TaylorF2 phase [72H75] and the 3PN amplitude (ZS [72] 
for compact binaries with comparable masses and spins 
(anti-)aligned with the orbital angular momentum. 

The FN expansion of the binding energy £ of such 
systems in the adiabatic approximation can be taken 
from the literature, see for instance |7TJ [TSHSU] and ref- 
erences therein. For the results shown here we include 
leading order and next-to-leading order spin-orbit effects 
[121 [SU [55] as well as spin-spin effects that appear at rel- 
ative 2PN order [8ll [83l [84] ; note that the square terms 
in the individual spins are valid only for black holes as 
discussed in 80, 83, 84j. The notation used in this section 
adopts unit total mass M = \ and G = c = 1. Each black 
hole is characterized by its mass nii and the magnitude 
of its spin 



\Xi\-m, 



1,2. 



(3.1) 



The spin vectors are (anti-)aligned with the orbital an- 
gular momentum L, where the sign oi L ■ Si defines the 
sign of Xi- With the aim of matching to available NR 
data, we use the FN spin definition that yields constant 
spin magnitudes [82l HSj . The quantity 



V 



mi TO2 
M2 



(3.2) 



is the symmetric mass ratio. The FN expansion is written 
in the dimensionless variable x, related to the orbital 



= ,,,2/3 



To 



angular frequency ui of the binary via x 

summarize the structure of this derivation, we start by 

giving the energy for the considered scenario as 



f = 



fe=0 



(3.3) 



where the coefficients e^ are listed in Eq. ( Al I 



The other ingredient needed to describe an inspiraling 
BBH as a sequence of quasi-circular orbits is the fiux J-", 
which we take at 3.5FN order including the same spin ef- 
fects as for the energy. We additionally take into account 
the 2.5FN correction of the flux due to the energy flow 
into the BHs, calculated in 1861. The flnal result is 



32 



J- = —TJ^x'^ J2 fk ■ 



,fc/2 



(3.4) 



fc=0 



where the coefficients fk are given in Eq. (A2). 



The energy loss of the system due to gravitational ra- 
diation is expressed as d£(t)/dt = —iFit), which trans- 
lates to an evolution equation for the orbital frequency, 
or equivalently 



dx •^{^) 

dt d£{x)/dx 



(3.5) 



Starting from (3.5), different waveform models can be 
constructed, for overviews see [501 [HZ]- For the purpose 
of the results shown here, we shall give the relevant ex- 
pressions in the frequency domain later and explicitly 
construct only the TaylorT4 approximant, which is ob- 
tained by expanding the right-hand side of Eq. (3.5) to 
3.5PN order 



dx 
'dt 



64 



-rjx 



k=0 



(3.6) 



with Ofc given in (A3). 



Note that the formal re-expansion of the denominator 
and the multiplication with the numerator in Eq. (3.5) 



also yields contributions to higher orders than those in 
Eq. (3.6). However, since 4PN and higher terms in flux 



and energy are not fully determined, the expressions one 
can compute for a^ with k > 7 are incomplete. The same 
applies to contributions of the spins at relative PN orders 
higher than 2.5PN. When we later use the TaylorT4 ex- 
pression (3.6) in this paper, we only expand it to 3.5PN 
order but keep all the spin terms that appear, i.e. incom- 
plete contributions in uq and ar are not neglected. Only 
if higher order spin corrections at 3 and 3.5PN order be- 
come available in energy and flux, the corresponding spin 
terms in the Taylor T4 (and TaylorF2) description can be 
completed. 

In order to construct an analytical formula of the wave 
signal in the Fourier domain, the stationary phase ap- 
proximation is commonly used to obtain the TaylorF2 
expression for the phase [72H75] . Below, we briefly reca- 
pitulate the steps towards the derivation of this approx- 
imant and provide the final result. 

Expanding the inverse of relation (3.5), dt/dx = 



— {d£ldx)/!F, allows for the analytical integration oit{x) 
The orbital phase (j) can be integrated via 



d<i) 
'dt 



3/2 
= LO = X ' 



d<i) 

dx 



= -x3/2 



d£{x)/dx 



(3.7) 



to obtain 4){x). This is the definition of the TaylorT2 
approximant. The {(., m) modes of the decomposition of 
the gravitational radiation in spherical harmonics can be 
approximated in the time domain by |76) 



/i«"(i)=A'^"(i)e 



-i7n(p{t) 



(3.8) 



and the transformation to the frequency domain is car- 
ried out in the framework of the stationary phase approx- 
imation 



/»'"(/) = 



/t^™(i) e'^'^'f'dt 



(3.9) 



^'"(i/)l 



27r 



i(f>{tf) 



^^r^u) 



(3.10) 



where i/ is defined as the moment of time when the in- 
stantaneous frequency coincides with the Fourier vari- 
able, i.e., mu}{tf) — 27r/. The phase in the frequency 
domain is given by 



r"'{f)^27rftf^mq)itf)- 



(3.11) 



Given t{x) and (f>{x) one can immediately change to the 
Fourier variable by 



X{tf) = [iu{tf) 



,2/3 _ /27r/ 



2/3 



(3.12) 



Starting from the energy ( 3.3 ) and flux ( 3.4 ) consequently 
leads to 

■ ^ (./)-5/3^a,(7r/)^-/3^ ^^-'^^ 



128r; 



k=0 



with the corresponding coefficients at of ( A4 1 . From 



(3.11) one realizes that, in fact, Eq. (3.13) is valid for 
all spherical harmonics with m = 2. The quantities io 
and 00 a-re arbitrary and arise as integration constants 
when calculating t{x) and (f>{x). When implementing this 
Fourier domain phase we also take into account the spin 
terms that appear after re-expanding at 3PN and 3.5PN 
order, although they are not complete. 

The time-domain amplitude of the gravitational wave 
was recently calculated at 3PN order by Blanchet et 
al. [7B]. We use the expression given by them for the 
^ = 2, TO = 2 mode in combination with the spin cor- 
rections provided in |26[ 177] . In our notation, the time- 
domain amplitude reads 



A'^ix) 



Si] X 



6 



5 ^^ 

k=0 



Aux'^l\ 



(3.14) 



where I?i, is the luminosity distance between source and 
observer and the coefficients Ak are given in ( A5 ) . 
From (3.10|) we see that, in order to construct the 



Fourier domain amplitude, an explicit expression for 
= (P'f^jdt^ = w is needed. In [77j this is done by re- 
expanding \/l/io using the same ingredients as those un- 
derlying the TaylorTn approximants. We may, however, 
look at w = (3/2) y^i: and choose one "preferred" pre- 
scription for X without re-expanding the quotient. Aim- 
ing at matching PN results to NR waveforms, we compare 
different possibilities of replacing x (namely by its Tay- 
lorTl and TaylorT4 description) and the re-expansion of 
the form 



OTT 

96^^ 



'11/4 



fe=0 



kX' 



fc/2 



(3.15) 
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FIG. 2. Different variants of constructing the PN Fourier 
amplitude in the stationary phase approximation for the 
equal-mass case. The labels explain how (Tr/cf))^'^ is treated 
in (3.101. The thick curve shows data obtained by a nu- 



merical simulation in full general relativity which begins at 
Af f ~ 0.008. The straight gray line illustrates the restricted 
PN amplitude, |ft^^|DL = 7ryV3(7r/)-'^/^ 



(see |77| ) with data of numerical simulations in full gen- 
eral relativity. The result in the equal-mass case can be 
seen in Fig. [2] Note that the transfer to the Fourier 
domain is completed by using (3.14) in (3.101 in combi- 
nation with (|3.12[). 



All variants of the 3PN Fourier amplitude agree rea- 
sonably well with the numerical relativity data roughly 
up to the frequency of the last stable circular orbit in the 
Schwarzschild limit, Mf = tt"! 6'^/^ w 0.022. Due to a 
comparable behavior even beyond this point we choose to 
construct the Fourier amplitude of our post-Newtonian 
model by using the TaylorT4 i (3.6). The same choice 
was employed e.g. in [88| . 



IV. MATCHING POST-NEWTONIAN AND 
NUMERICAL RELATIVITY WAVEFORMS 

A. Basic notions 

The basic criteria for evaluating the goodness of fit for 
the hybrid waveform require a notion of distance between 
two GW signals h{t) and h'{t). The simplest notion is 
the distance in the least-squares sense over an interval 
ti<t<t2 [37H39] 



Jtut2 



{h,h') 



\h{t)-h'{t)f dt. 



(4.1) 



This can be used for the numerical relativity /iNR(i) and 
the post-Newtonian waveform /ipN(t), with the interval 
[ii,i2] being chosen so that both waveforms are reason- 
ably good approximations (in a sense to be quantified 
later). Thus, the PN waveform is taken up to ^2 and the 



NR waveform is taken to start at ti, and they overlap 
within the interval [ti , t2] . 

Let us consider the frequency domain equivalent. Our 
convention for the Fourier transform of a signal x{t) is 



(/)= / x{t)e'-^f'dt. 



(4.2) 



One needs to be careful in converting the time interval 
[^1,^2] to a frequency interval [/i,/2]- In principle, the 
Fourier transform is "global" in time; signals that have 
compact support in time cannot have compact support in 
frequency, and vice versa. However, for the binary black 
hole waveforms (prior to the ring-down stage) that we 
are considering, the frequency always increases in time, 
so that we can sensibly associate a frequency interval 
[/ij /2] with a given time interval [ti, ^2]- For these wave- 
forms, we can consider the above distance definition in 
the frequency domain: 



/2 



SfuhiKh')^ h{f)~h'{f) 



/i 



df. 



(4.3) 



We shall use such a norm (applied to the phase) for con- 
structing the hybrid waveform. 

When evaluating the goodness of a hybrid waveform 
for a particular detector, we need to consider detector- 
specific inner products, which are convenient to describe 
in the frequency domain. Let Snif) be the single-sided 
power spectral density of the noise in a GW detector 
defined as 

E[h{f)h*if)]^lsM)S{f-f). (4.4) 

Here n{t) is the detector noise time series with h{f) its 
Fourier transform, * denotes complex conjugation and 
E refers to the expectation value over an ensemble of 
independent realizations of the noise, which is assumed 
to be a zero-mean, stationary, stochastic process. This 
equation implies that data at different frequencies are 
independent, and is one of the reasons why working in the 
frequency domain is so useful in data analysis. The time 
domain description of the noise is more complicated; n{t) 
and n{t + T) are in general not independent; E,[n{t)n{t + 
t)] is generally non-zero. For stationary noise this is a 
function C{t) only of r, and is related to S'„(/) via a 
Fourier transform (see e.g. [521 )■ 

Given S'„(/), we use the following definition of an inner 
product between two signals x{t) and y{t) 



{x\y) = 4Re 



Snif) 



df, 



(4.5) 



where x{f),y{f) are the Fourier transforms of x(t),y{t) 
respectively. This inner product is appropriate for Gaus- 
sian noise and forms the basis for matched filtering (see 
e.g. [9^). It can be used to define a suitable notion of 
distance between two signals h{t) and h'{t) as {Sh\Shy/^, 
where Sh{t) = h'{t)-h{t). 



The distinguishability between h{t) and h'(t) in the 
presence of noise can be understood with the following 
construction. Following [HI], we define a 1-parameter 
family of waveforms which interpolates linearly between 
h[t) and h'{t) as 



h"{t;X) ^h{t) + Xdh{t). 



(4.6) 



We obviously have h"{t;0) = h{t) and h"{t;l) = h'{t). 
The question of distinguishability between h{t) and h'{t) 
now becomes one of estimating the value of A [for the 
extended signal model h"{t;X)] in the presence of noise. 
If we use an unbiased estimator for A, the variance a^ of 
the estimator is bounded from below by the Cramer-Rao 
bound (see e.g. [Hj) 



al > {dh\Shy 



(4.7) 



This can be a useful bound for large SNRs, which is in 
fact what we are interested in here; it is easier to dis- 
tinguish between two loud waveforms and demands on 
the waveform model are correspondingly more stringent. 
Thus, a useful condition for being able to distinguish 
between the two waveforms is ca < 1. If h(t) is the 
true waveform and h'{t) our approximation to it, then 
we say that h'{t) is a sufficiently accurate approximation 
if {6h\Sh) < 1. Assuming that (h\h) ~ (h'\h'), it is clear 
that {Sh\Sh) ex p^ where p = {h\hy/^ is the optimal SNR. 
Hence, as we just remarked, the two signals are easier to 
distinguish when the detector is more sensitive, or when 
the signal amplitude is larger. It will be convenient to 
normalize the norm of 6h and write this distinguishability 
criterion as 



XiSh\Sh) > \ 



(4.8) 



Thus, for a given detector, we choose a reasonable guess 
po for the largest expected SNR and we compute the nor- 
malized distance between the two waveforms {Sh\Sh) / p^ . 
If this exceeds 1/pq, then we consider that the detector 
is able to distinguish between the two waveforms. 

If we are interested in the less stringent requirement of 
detection rather than in strict distinguishability, then a 
sufficient condition is El] 



1 



iSh\Sh) < 2e, 



(4.9) 



where e is the maximum tolerated fractional loss in SNR. 
More explicitly: if h{t) is the exact waveform and h'{t) 
an approximation thereof, then the approximation is po- 



tentially useful for detection purposes if (4.9) is satisfied 



for an appropriate choice of e. If we are willing to accept 
e.g. a 10% loss in detection rate, then a suitable choice 
is e « 0.10/3 « 0.03 (corresponding to sources uniformly 
distributed in space). This value of e does not take into 
account the additional loss in SNR due to discrete tem- 
plate banks used in realistic searches. In practice one 
might need to choose e an order of magnitude smaller 



than this so that the total fractional loss in SNR remains 
acceptable. Since a more precise value is pipeline de- 
pendent, we shall ignore this caveat and use e = 0.03 
as a convenient reference; the reader can easily scale the 
results of this paper appropriately for different choices. 

A useful way to describe the efficacy of approximate 
waveform models is through the concepts of effectualness 
and faithfulness introduced in [HS]. Let h\{t) be the ex- 
act waveform with parameters A and the approximate 
waveform model be h^^(t). The ambiguity function is 
defined as the normalized inner product maximized over 
extrinsic parameters 



A{X, A') = max 



(hxlhT) 



*o.*o vw^ihwnwy 



(4.10) 



where to is the time offset between the two waveforms, 
and (/)o is the initial phase. Performing a further maxi- 
mization over the parameters A' of the model waveforms, 
we define -4(A) = max^/ ^(A, A'). If A{X) exceeds a cho- 
sen threshold, e.g. 0.97, then the waveform model /i'^pp 
is said to be effectual. Effectual models are sufficient for 
detection. In order to be able to estimate parameters 
we also need the model to be faithful. This means that 
the value of A' which maximizes .4(A, A') should not be 
biased too far away from A. 



B. Issues in matching PN with NR 

It is useful at this stage to discuss some of the issues 
that arise in combining PN and NR results. The dis- 
cussion here will be short and incomplete, and the topic 
merits an in-depth investigation that is beyond the scope 
of the present work. Our immediate aim is simply to spell 
out some of the reasons why black hole parameters in the 
PN and NR frameworks may not necessarily refer to the 
same physical quantities. One should therefore not be 
surprised that when combining NR and PN waveforms, 
it might become necessary to vary the intrinsic black hole 
parameters as well. This is not to say that cither PN or 
NR use incorrect definitions for black hole parameters, 
both frameworks are in fact consistent within their do- 
mains of applicability. The point rather is that the two 
formalisms are quite different when viewed as approxima- 
tion schemes to general relativity, and these differences 
might need to be taken into account depending on the 
accuracy requirements for the matching. 

Since PN and NR are both used to address the BBH 
problem, one could imagine starting with the two black 
holes very far apart, evolve them using appropriate PN 
equations of motion and compute the resulting wave- 
forms. As one gets close to the merger, terminate the 
PN evolution and use this end-point to construct initial 
data for the full NR simulation which then evolves the 
black holes through the merger and ringdown. However, 
the formalisms and methods employed in the two cases 
are radically different and there are potential difficulties 
in carrying out this procedure. 
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The PN formalism is based on a perturbative expan- 
sion in powers of the small parameter e = w/c, where 
V is the orbital velocity and c is the speed of light. In 
the usual formulations, PN theory uses a point-particle 
description of the black holes, and their parameters can 
be viewed as effective parameters which couple in the ap- 
propriate manner with the external background gravita- 
tional field (see e.g. [MlinSj)- The goal of PN theory is to 
find a 1-parameter sequence of solutions to the field equa- 
tions gf^^ to any specified order in e. It has recently been 
shown rigorously (HS] that, in the cosmological setting 
with gravitating perfect fluids, the 1-parameter family of 
solutions exists and admits an expansion in e" to any or- 
der. While similar results in the asymptotically-flat case 
are not yet available, it is certainly reassuring to know 
that PN works well in this non-trivial setting (in fact, it 
can be persuasively argued that the cosmological setting 
is more relevant to GW observations than strict asymp- 
totic flatness). The errors in PN waveforms are then due 
to the systematic differences between the true waveform 
and the asymptotic series expansion in e" truncated at a 
flnite order, and this error depends on which particular 
PN expansion one chooses to use. 

In contrast, numerical relativity is based on the 3-f 1 
formulation of general relativity as an initial value prob- 
lem, and one solves the resulting partial differential equa- 
tions numerically. The GW signal is typically measured 
at a large, but finite, coordinate distance from the source, 
and encoded in ^4, the frame-dependent outgoing trans- 
verse component of the Weyl tensor component. The 
data from multiple coordinate radii are extrapolated to 
asymptotic distances, or evaluated at null infinity in the 
case that characteristic extraction is used tSHl EH] • For 
a given physical configuration (choice of masses, spins, 
separation etc.), one specifies the initial data consisting 
of the spatial metric and extrinsic curvature of the ini- 
tial spatial slice. The physical initial data parameters 
should be chosen to be as compatible as possible with 
the spacetime computed in the PN formalism, and sig- 
nificant progress has been made in this regard [23 HH] • 

The black holes here are not point particles but rather 
black hole horizons. The parameters of the black hole 
are often computed as integrals over the apparent hori- 
zon, and in most cases the parameters used in construct- 
ing the initial data are also useful approximations to the 
true ones. There are, however, possible systematic er- 
rors. For example, if we are using the quasi-local horizon 
definitions, an important requirement is that the hori- 
zon should locally be approximately axisymmetric. The 
methods for finding the approximate symmetry vectors 
have become increasingly accurate and reliable [55Ul03j . 
However, it should be kept in mind that the assumption 
of approximate axisymmetry is expected to become in- 
creasingly worse closer to the merger. Furthermore, the 
very use of apparent horizons is gauge dependent; using a 
different time coordinate will lead to a different set of ap- 
parent horizons and possibly also different values of the 
parameters. In the inspiral phase when the horizons are 



sufficiently isolated this gauge issue is not expected to be 
a problem, but as we get closer to the merger (this has 
not yet been quantified), the variation in the parameters 
due to gauge choices could become significant |104j . 

Let us elaborate a little more on the spin. Most post- 
Newtonian treatments are based on the equations of mo- 
tion derived in [1051 1106] . The starting point is the spin 
tensor S''"' constructed from moments of the stress en- 
ergy tensor T'"^. Since 5'"' has potentially 6 non-zero 
independent components, the system for the 4 equa- 
tions of motion V^T'^" = is over-determined. One 
thus imposes the additional spin supplementary condi- 
tions such as S^^^pv = or S^^u^ = with p^ being 
the 4-momentum and Ui, the 4-velocity. These different 
conditions lead to physically different equations of mo- 
tion and trajectories [107] . On the other hand, for black 
holes in NR, a common method for evaluating spin em- 



ploys the formalism of quasi-local horizons [108] . The 
final result for the magnitude of the horizon angular mo- 
mentum is an integral over the apparent horizon S: 



J =-- 



1 

8^ 



Kf,^<j)^'dS'' , 



(4.11) 



where K^^^ is the extrinsic curvature of the Cauchy slice, 
(j)^ is a suitable approximate axial symmetry vector on 
S [MHlOSj . and dS^ is the area element on the appar- 
ent horizon. The direction of the spin is harder to find, 
but some approximate methods are available [lOOL I109| . 
There is yet no detailed study of possible analogs of the 
spin supplementary conditions in this formalism, or on 
the equations of motion for horizons with a given set of 
multipole moments. For a horizon with area A and spin 
magnitude J, the mass is given by the Christodoulou 
formula 



TO 



A 
16^ 



47rJ2 
A 



(4.12) 



Hence, uncertainties in spin can also lead to uncertainties 
in the mass. 

As long as we are dealing with just the numerical or PN 
waveforms by themselves, these small effects in the defini- 
tions of mass and spin are not important for most appli- 
cations. In fact, we can treat them as just convenient pa- 
rameterizations of the waveform without worrying about 
their detailed physical interpretation. However, when we 
wish to compare the results from frameworks as differ- 
ent as PN and NR this may no longer work. Depending 
on the details of the matching procedure, systematic dif- 
ferences between the various definitions might need to 
be taken into account, or at the very least they should 
be quantified. If a particular case requires matching a 
very long PN portion (depending on the total mass and 
the lower- frequency cutoff of a particular detector), then 
even a small change in the PN parameters at the match- 
ing frequency can translate into a large phase difference 
at lower frequencies. One valid approach is to not as- 
sume a priori that the PN and NR parameters are equal 
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to each other but rather, for a given numerical waveform, 
we search over PN waveforms in a particular PN approx- 
imant and find the best fit values. Finally, given that the 
NR and PN parts have different values of physical pa- 
rameters, it is a matter of convention what values are to 
be assigned to the hybrid. The values of the parameters 
in the early inspiral are a convenient and astrophysically 
relevant choice. 



C. An illustration for non-spinning systems 

Let us now move to a concrete case of constructing 
hybrid waveforms, considering the non-spinning Llaina 
waveforms, i.e. data set #7 in Table [ij Recall that this 
data set consists of two waveforms with non-spinning 
black holes with mass ratios 1:1 (used in left and central 
panels of Fig. [6| and 1:2 (Figs, [s] |4J |5] and right panel 
of IgI) . Since these waveforms are calculated using the 
Llama code with extraction at future null-infinity with 
the Cauchy-characteristic method for the equal mass 
case, or well into the wave-zone for the 1:2 case, we are 
confident that systematic effects of waveform extraction 
are small. Even for these waveforms, based on the discus- 
sion above, in principle we should not rule out a small 
mismatch in the values of the spin (and perhaps also 
eccentricity) between the NR and PN waveforms. For 
simplicity, let us consider only the possibility that the 
symmetric mass ratio t] could be different, and restrict 
ourselves to non-spinning black holes and zero eccentric- 
ity. We would like to match the Llama waveforms with 
the frequency domain PN waveforms discussed in Sec. |III| 
with the values of the spins set to zero. The total mass 
M sets the scale for the time (and frequency) ; in addition 
we have the extrinsic parameters for the time offset and 
initial phase to ^nd (J)q. Furthermore, we only consider 
the £ — m — 2 mode, so that the PN waveform is of the 
form h^^ {AIf;(f>o,to,rip^) in the frequency domain. 



1. Fitting errors 

For a given NR waveform hjq^{t) we consider a time 
window {to-, to + Ai) or, alternatively, in the frequency 
domain the matching region consists of a lower starting 
frequency /^ and a width A/. We match the two wave- 
forms in a least-squares sense by minimizing the phase 
difference in Fourier space 



pfL+Af 

S — min / 

to,0O,'7PN J fj^ 



LogioA<^o 



l<5(?!'(/;?7NR,??PN,io>o)| Mdf , 
S<t){f) = (/'nr(/; ?/nr) - 0pn(/; to, (po, i7pn) ■ (4.13) 



We optimize S over all allowed time and phase shifts, i.e. 
{to, 00 ) J and the PN intrinsic parameters ApN- Given the 
previous discussion on the possible differences between 
the intrinsic parameters A in the PN and NR frameworks. 




0.009 



0.012 0.015 O.OIS 

Mf/, fitting window 



FIG. 3. A contour plot for the fitting error A(j)o in the 
(/l,A/) plane. Here rj is kept fixed to the NR value and 
we optimize over 00 and to- 



here we have distinguished between the intrinsic param- 
eter rj (3.2) appearing in /ipN and /inr- Note that we are 



not only neglecting spins and eccentricity but also assume 
A/pN = Mnr = M. Future analyses should successively 
drop these simplifications. 

Let us now consider the choice of the optimal match- 
ing window {JltIl + A/), and the best fit values of 
{4>o,to,r]pf^). For each window, the least squares pro- 
cedure gives a best fit value ?7pn = v{fL,^f) and l-a 
error estimates Ary, A(j)o, Ato. Our principle for choosing 
{fL,Af) is to pick the one for which the quality of fit 
between the NR and PN waveforms is the best, i.e. to 
minimize the fitting errors. 

We first fix r;pN = t^nr, choosing the 1:2 waveform, and 
consider fitting for (0Oi^o)- The result for A0o is shown 
in Fig. [3] as a contour plot in the (/l. A/) plane. There 
are clearly multiple best-fit islands but we already see 
that the optimal window choice turns out to be a long 
frequency width starting at low frequencies, or a rela- 
tively short window starting closer to the merger. Re- 
garding the increasing error PN most likely introduces 
towards higher frequencies, we prefer using an early and 
long matching window. Though we do not show it here, 
the result is similar for the time offset to- 
ll is more interesting instead to generalize this and 
allow all three parameters (?7pn, 0Oj ^0) to vary. The main 
result is displayed in Fig. [4j which shows contour plots 
of the fitting errors A77, A0o and Atp in the (/l. A/) 
plane. There are now clear and consistent minima for all 
errors and thus a clear best choice for Jl and A/. At 
this optimal choice, we see that we can fit r], (po and to to 
better than 10^'^, 0.06 and 0.15M, respectively. Apart 
from the error Arj, the actual best fit value rj is also of 
great interest. Fig. [5] shows the value of 77 as a function of 
the start frequency of the matching window /^ and A/. 
The X-axis on this plot is the start point of the fitting 
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FIG. 4. Dependence of the fitting errors in rj, 4>o and io on the frequency window {fL,Af). Note that there is a clear 
choice of (/l, A/) ~ (0.0093,0.014) that optimizes the fit between the NR waveform and the PN waveforms with different 
rj. For a binary of total mass 10 Mq, for which the last stable orbit happens at 440 Hz, this corresponds to frequencies 
(/l, A/)|io A/g ~ (189, 284) Hz. This indicates that the optimal window for matching should start at the lowest reliable 
frequency available from the NR waveform and extend roughly up to the last stable orbit, which usually quantifies the point 
when the PN approximation starts to break down. Moreover, the plot on the left shows that the accuracy in i] decreases 
slowly with different choices of (/l. A/), assuring that small changes in these values do not lead to large errors in the hybrid 
construction. At the best fit point, the accuracy in r; by this fitting procedure is better than 10"'^. 
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FIG. 5. Best fit value of ry as a function of the start fre- 
quency /l of the matching window for the waveform which 
corresponds nominally to a mass ratio 1:2, i.e. t^nr ~ 2/9 = 
0.222 . . .; this is shown by a horizontal dashed line. The ver- 
tical dashed line at M/_l = 0.009 is the start frequency of the 
NR waveform. A rectangle highlights the region of minimal 
fitting errors from Fig. [4] We see that the best determined 
values of rj are clearly less than t/nr- 



window /l, and the color bar indicates A/. The most 
trustworthy values correspond to the optimal choice of 
(/l, A/) obtained in Fig. |4J we indicate the union of all 
three minimal-error islands as a rectangle in Fig. [SJ 

To summarize, from Figs. |4] and [5] we deduce that, if 
we were to ignore ?7nr (the value that the numerical sim- 
ulation nominally assumes) and simply try to find the 
best fit with the PN waveforms described in Sec. |III[ 
then we can clearly estimate the best matching region 
(/l, /l + A/) and a best fit value yypN = rj ± Ai]. This 



procedure illustrates a trade-off between trying to match 
at early frequencies, where our PN model is more reliable 
and having a sufficiently long fitting window, in which a 
considerable frequency evolution leads to an accurate es- 
timate of the fitting parameters. The difference between 
?7nr and rypN for this case is seen to be ~ 10%. This 
by itself does not say that the uncertainty in rj is 10% 
because as we shall soon see, the uncertainties in the hy- 
brid waveform are dominated by the uncertainties in the 
PN model. In other words, the NR waveform is closer to 
the true physical waveform within the matching window, 
and we should not actually use the best fit value of ?7pn 
to construct the hybrid. 



2. Accuracy of the hybrid waveform 

Later we shall show a phenomcnological fit for the hy- 
brid waveform and we shall claim that the fit reproduces 
the hybrid waveform sufficiently accurately. Here we first 
ask whether the hybrid waveform is itself sufficiently ac- 
curate subject to various errors. The basic criteria for 
evaluating this is the notion of a distance between two 
signals whose difference is 6h, as given in Eq. (4.8). For 
two signals h and h' , we shall consider the normalized 
distance squared {6h\Sh)/p'^, where p is calculated from 



our best model (3PN amplitude, 3.5PN TaylorF2 phase 
combined with highest resolution NR waveform). Now 
the total mass M becomes important. Previously, when 
we looked at the least square fits in Eq. (4.13), the to- 
tal mass appeared just as a scale factor. However, in 
the inner product Eq. (4.51, the power spectral density 



Sn{f) sets a frequency scale, and the value for {Sh\Sh) 
becomes mass-dependent. We shall consider two design 
noise curves. Initial and Advanced LIGO |110lllll] . We 
are then addressing the question of how different our hy- 
brids would be if we were to use a slightly different result 
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on either the NR or PN side. 

On the NR side, we first consider data computed at 
different resolutions. The Llemia waveforms for the equal- 
mass case have been computed at low, medium and high 
resolutions corresponding to spacing h = 0.96, 0.80 and 
0.64 on the wave extraction grid. The finest grid, i.e. 
the grid covering the black hole, has a resolution of 0.02 
for the finest resolution. This is scaled by 0.80/0.64 and 
0.96/0.64 for the medium and low resolution runs respec- 
tively. We combine these waveforms with the TaylorF2 
model from Sec. |III| by using the optimal matching win- 
dow discussed around Fig. |3] and ?7pn = ?7nr- The result 
is shown in the left panel of Fig. [6J Hybrids constructed 
with medium- and high-resolution waveforms would be 
indistinguishable even with Advanced LIGO at a SNR 
of 80 over the considered mass range. Thus, we conclude 
that the numerical errors related to a finite resolution are 
not relevant in the hybrid construction process. 

The uncertainties increase when comparing NR data 
produced by different codes. Similar to the analysis of 
different resolutions we calculate the distance of hybrid 
waveforms for non-spinning black holes with mass ra- 
tio 1:1 and 1:2. Results from data set ^^1 and ^^^8 (see 
Table IT]) were used, and the distance plot in the cen- 
tral panel of Fig. [6] shows that the 1:2 waveform would 
be distinguishable for Advanced LIGO at SNR 20 for 
a total masses between ~ SOMq and ~ 65Mq. Note 
that these errors are dominated by our matching to PN 
which possibly yields different fit parameters for the PN 
model and therefore amplifies small differences in the NR 
data. Towards higher masses, the influence of this match- 
ing decreases as well as the distance of both waveform. 
However, as we shall show next, all these errors are still 
small compared to the intrinsic uncertainties introduced 
by PN and they do not matter for Initial LIGO. If we care 
only about detection with a minimal match e — 0.03 [see 
Eq. (4.9)], we have even less to worry about. 



The errors on the PN side turn out to be much more 
important. The right panel of Fig. [6] illustrates the effect 
of using different PN approximants combined with the 
same SpEC equal mass simulation. We first use the fit- 
ting window discussed above, although the exceptionally 
long SpEC waveform would allow a much earlier match- 
ing. The dashed curve shows the difference in the hybrid 
waveforms when we match the 3PN or 3.5PN phase fol- 
lowing the TaylorF2 frequency domain approximants de- 
scribed in Sec. HI (the amplitude is taken at 3PN order 
in both cases). We see that the difference between these 
hybrids becomes significant even for Initial LIGO at SNR 
of 8 between a total mass of ~ 5Mq and ~ 35Af0 . Sim- 
ilarly, the differences between the F2 and Taylor Tl & 
T4 approximants are also significant. For detection with 
e = 0.03 [see Eq. ( |4.9[ )], we need to look at the horizontal 
line with {Sh\6h)/p"' — 0.06 in the right panel of Fig. 6] 
Both in the 3PN/3.5PN distance (dashed line) and the 
TaylorTl/TaylorF2 comparison (upper black solid line), 
there is a small range of masses for which the difference 
between the hybrids would matter even for detection. 



As a reference, we make use of the fact that the numer- 
ical data #9 (Tablell]) actually contains physical informa- 
tion to frequencies considerably smaller than the match- 
ing window used for our hybrid production. We therefore 
match the TaylorF2 phase at 3PN and 3.5PN order also 
with a much earlier fitting window (roughly a factor of 
2 lower in frequency). The right panel of Fig. p] shows 
that the difference indicated as "early match" remains 
undetectable for a larger range of total masses. Expand- 
ing such studies may be used to quantify the necessary 
length of NR waveforms and estimate to which frequency 
standard PN results can be used in hybrid waveform con- 
structions. 

Having carried out this study of errors for non-spinning 
waveforms, we can now draw some conclusions for the 
aligned-spin case. In principle, the procedure outlined 
here remains valid; we should search over not only 
{ri,to,(j)o}, but now also over the spins {xi,X2}- We 
would not expect the results to be better than shown here 
for non-spinning waveforms because (i) we are adding 
two more parameters and (ii) the waveforms #1-4 are 
expected to have more wave-extraction systematic errors 
than the Llama results considered here. Most impor- 
tantly, as we have just seen, the intrinsic errors in PN 
are more significant whereas the numerical accuracy is 
not the bottleneck. The intrinsic parameter biases in PN 
also show up when different PN models are compared 
with each other. An extensive comparison of different PN 
models is made in [87"; this paper quantifies the mutual 
effectualness and faithfulness of the different PN mod- 
els and shows that errors of ^ 20% are not uncommon 
for Advanced LIGO. The less than 10% discrepancy in rj 
shown in Fig. [s] is thus entirely consistent with the differ- 
ences between different PN models. To address this, one 
needs either improved PN models or a greater variety of 
longer NR waveforms such as the long SpEC simulation. 

As a simplification, in what follows below we will 
choose the matching window based on maximizing over 
the extrinsic parameters (io,0o) motivated by Fig. [3J In 
that figure, we observe the best fit region extending di- 
agonally from MAf w 0.013 on the y-axis, to the bot- 
tom right corner. It turns out that for this diagonal, 
the upper frequency of the window does not vary much, 
0.020 < MJl + MAf < 0.024, and we shall use this 
fact below for constructing hybrid waveforms for aligned 
spinning systems. 



D. Construction of hybrid waveforms for aligned 
spinning systems 

Let us now proceed to the construction of a hybrid 
waveform model for non-precessing, spinning systems 
with comparable mass. Again, the waveforms described 



in Sec. HI will be the basis for our model at low frequen- 
cies corresponding to the inspiral stage. On the other 
hand, the NR simulations described as data-sets #1~3 
in Table IT] contain physical information for frequencies 
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FIG. 6. Distinguishability of hybrid waveforms that have been constructed varying some of the hybrid ingredients at a time. 
When indicated, the black/grey color code denotes that Initial/ Advanced LIGO design curves have been used for the distance 
calculation. The horizontal lines are the lines of constant SNR (in fact it is 1/SNR'^); if the distance measure goes above them, 
then the waveforms can be distinguished from each other. The left panel shows the effect of constructing hybrids from Llama 
equal-mass waveforms at different resolutions. We consider the difference between the high-medium resolution waveforms, and 
the high-low waveform resolutions. The central panel shows the effect of using NR waveforms produced with either BAM or 
Llama codes. The solid lines indicate the normalized distance in the equal-mass case, dashed lines show the case of mass-ratio 
1:2. The highest available resolution was always used. The panel on the right displays Initial LIGO's ability to distinguish 
hybrid waveforms constructed from different PN approximants. This plot shows that the hybrids are not sufficient for detection 
at the e = 0.03 level [Eq. (4.9l] only for a small range of masses. "Early match" is a reference for matching 3PN or 3.5PN F2 



at early frequencies to the long equal mass SpEC waveform. 



above Mf « 0.009. We will refer to Fig. [3] to justify our 
choice of an overlapping window at M f E (0.01, 0.02). 

Once this interval is fixed, we now carry out the follow- 
ing matching procedure for all NR simulations of data- 
sets #1-3: PN and NR phases are aligned by fitting the 
free parameters tg and cjjq in Eq. (3.13); with a standard 



root-finding algorithm (starting at the mid point of the 
fitting interval) we find a frequency /$ where PN and 
NR phase coincide and construct the hybrid phase con- 
sisting of TaylorF2 at / < /$ and NR data at / > /$. 
An analogous procedure is applied to the amplitude, but 
in this case there is no freedom for adjusting any param- 
eters. Hence, we use an educated guess for the matching 
frequency (compatible with that for the phase) and find 
the root /^ where the difference of PN and NR ampli- 
tude vanishes. The hybrid amplitude consists of PN data 
before and NR data after Ja- Small wiggles in the NR 
amplitude, due to the Fourier transform, do not affect 
the phenomenological fit significantly. The most impor- 
tant ingredient for arriving at an effectual model is the 
phase. 

Fig. [7] illustrates the above- described hybrid construc- 
tion method for matching PN and NR data in the fre- 
quency domain. The procedure does not require any re- 
sizing the PN or NR data and allows for the construction 
of waveforms containing all the information from the Tay- 
lorF2 approximant at low frequencies and input from the 
NR simulations for the late inspiral, merger and ring- 
down. The resulting hybrid PN-NR data cover a part 
of the parameter space corresponding to equal-valued, 
(anti-) aligned spins for 0.16 < rj < 0.25 and constitute 
the "target" waveforms to be fitted by the analytical phe- 
nomenological model described in Section Iv] 



V. PHENOMENOLOGICAL MODEL 

In this section we present the phenomenological model 
developed in order to fit the hybrid PN-NR waveforms of 



Section IV to an analytical formula. A geometric descrip- 
tion of the procedure for constructing phenomenological 
waveforms parametrized by just the physical parameters 
is detailed in [39, . and here we just summarize it. Let 
Ai be the space of intrinsic physical parameters that we 
are interested in. In the present case, this is the four- 
dimensional space of the component masses and spins 
A ~ {M, T], Xi, X2}- For each point A in M, let h{t; A) be 
the true physical waveform that we wish to approximate; 
in particular we consider only the dominant £ — m = 2 
mode in this paper. Furthermore, as in [40] . we model 
the spin effects using a single parameter x defined as 



1 



1 



X 



-Xi 



-X2, 



(5.1) 



where 6 = {mi — 1112) /M. As mentioned in Sec. pj this is 
justified from PN treatments of the inspiral 53 and from 
numerical simulations of the merger [54j which considers 
equal mass systems. In these works, it is found that 
the dominant spin effect on the waveform is from the 
mass-weighted total spin of the system. On the PN side, 
this can be further justified by looking at the expressions 
for the PN phase and amplitudes given in the appendix 
[Eqs. ( A4|A5 )]. Consider the phase and amplitude terms 
as polynomials in 77 and retain only the 0{r]^) terms. 
At this lowest order, the amplitude and phase are seen 
to depend only on x- We can therefore hope that this 
single parameter captures the main effects of the black 
hole spins, at least for the purpose of constructing an 
effectual model. In fact, whenever we incorporate pure 
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FIG. 7. Illustration of the method for constructing PN-NR hybrid waveforms in the frequency domain. The data corresponds 
to an equal-mass binary with aligned spins Xi — Xa = —0.25. The left panel shows the amplitude and the right panel displays 
the phase of the dominant (. = 1,m = 2 mode of the GW complex strain /i(/). The green dotted lines correspond to the 
TaylorF2 PN approximant and the red dot-dashed curve is the NR data. The hybrid waveform is depicted in solid black and 
the matching points for amplitude and phase are indicated with a dashed line. 



PN contributions in our final model, we use them with 
X\—Xi— X- It is however important to note that this is 
only an approximation; while it sufSces for our purposes 
(i.e., in constructing an effectual model) it will need to 
be refined as more faithful models are required. The 
degeneracy in the space of aligned spins will be studied 
in greater detail in a forthcoming paper [112] , but here we 
look to construct a phenomenological model using only 
(M, 77, x) as the physical parameters. 

We start with some known signals in this parameter 
space at N points Ai, A2, . . . , Aat. We take these known 
signals to be the hybrid waveforms whose construction 
we described earlier. Here the NR waveforms are the 
BAM waveforms of data sets #1-3 summarized in Table |l) 
and the PN model is the 3.5PN frequency domain model 
for aligned spins described in Sec. |III| Given the finite 
set of hybrid waveforms constructed from these ingre- 
dients, we wish to propose a phenomenological model 
^phcn(i;A) that interpolates between the hybrid wave- 
forms with sufficient accuracy. In constructing this phe- 
nomenological model, it is convenient to work not with 
the physical parameters A, but rather with a larger set of 
phenomenological parameters A, which we shall shortly 
describe. If A^ is the space of phenomenological pa- 
rameters, then we need to find a one-to-one mapping 
M. ^ M denoted A (A), and thus the subspace of M. cor- 
responding to the physical parameters. As the end result 
of this construction, for every physical parameter A, we 
will know the corresponding phenomenological parame- 
ter A(A) and thus the corresponding phenomenological 
waveform /iphcn(^; A(A)). 

Following the construction procedure of Section |IVD[ 
we split our waveforms in amplitude and phase, both of 
which shall be fitted to a phenomenological model 



For both the amplitude and the phase of the dominant 
mode of the GW radiation, we make use of the insights 
from PN and perturbation theory for the description of 
the inspiral and ringdown of the BBH coalescence, re- 
spectively, and introduce a phenomenological model to 
complete the description of the waveforms in the merger. 



A. Phase model 

The PN approach for the GW radiation based on the 
stationary phase approximation, introduced in Eq. (3.13) 
of Section 



III| (used with ig = 0o = and xi = X2 = x)j 
gives an adequate representation of the phase of the dom- 
inant mode during the adiabatic inspiral stage of the 
BBH coalescence V'spaI/)- ^^ the system transitions to- 
wards the merger phase, it is expected that further terms 
in the expansion are required to capture the features of 
the evolution. With this ansatz in mind, we propose a 
pre-merger phase V'pm(/) of the form 






"3/ 



-1/3 



+04 + 05/2/3 + ^6/), (5.3) 



where the a^ coefficients are inspired by the SPA phase, 
redefined and phenomenologically fitted to agree with the 
hybrid waveforms in the region between the frequencies 
0.1/rd and /rd, which depend on the spins and masses of 
the black holes in the form explained below in Eq. (5.5 1. 



Ve„(/)= Ven(/)e'*^^""(^). 



(5.2) 



Note that 0.1/rd roughly corresponds to the starting fre- 
quency of our NR simulations. 

As for the post-merger phase, the Teukolsky equa- 
tion |113j describes the ringdown of a slightly distorted 
spinning black hole. The metric perturbation for the fun- 
damental mode at large distances can be expressed as an 
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FIG. 8. Fitting procedure for the amplitude, applied to the 
equal-mass, non-spinning case. The 71 term of Eq. ( |5.11[ ) 
is introduced to follow the behavior of the amplitude in the 
pre-merger regime whereas the Lorentzian curve correctly de- 
scribes the post-merger. The two pieces are glued together in 
a smooth manner using tanh-windows. 



exponential damped sinusoidal 



/»?L(0 = 



■Aring-'" ^_7r 



f-R-ot/Q p-2-Kif-R_T> t 



(5.4) 



where M is the mass of the ringing black hole, D^ the 
distance from the source, and Q and /rd correspond to 
the quality factor of the ringing down and the central fre- 
quency of the quasi-normal mode. These can be approx- 
imated with an error < 2.5% in the range a € [0, 0.99] by 
the following fit [114] 



/RD(a,M) 
0(a) 



1 



2ttGM 

qi +92(1 



[fci + fc2(l-a)''^ 



(5.5) 
(5.6) 



where h = {1.5251,-1.1568,0.1292} and q, = {0.7000, 
1.4187, -0.4990}, i ^ 1, 2, 3 as given in Table VIII of [TT4] 
for the {I, m, n) — (220) mode. The review [115] presents 
a ftrll description of quasi-normal modes. The quantity 
aM'^ is the spin magnitude of the final black hole after 
the binary has merged, which can be inferred from the 
spins of the two black holes. In our case, we use the fit 
presented in |57| . which maps the mass-ratio and spins 
of the binary to the total spin a of the final black hole. 



The analytical treatment of the ringdown (5.4) moti- 



vates a linear ansatz for the post-merger phase '0rd(/) 
of the form 



V'rd(/) = /3i + /32/. 



(5.7) 



The /3i 2 parameters are not fitted, but obtained from the 
pre-merger ansatz (5.3 1 by taking the value and slope of 
the phase at the transition point /rd- The transition 
between the different regimes is smoothened by means of 
tanh-window functions 



'fo 



1 ± tanh 



4(/ - /o) 



(5.8) 



to produce the final phenomenological phase 

^phon(/) = V'ipAWj^ + V'pmW/i wJ^ + V'RD^/a : (5-9) 

with /i = 0.1/rd, /2 = /rd; here we have used d = 0.005 
in the window functions w^ . Roughly, these two tran- 
sition points respectively signal the frequencies at which 
our NR simulations start and the point at which the bi- 
nary merges, and have been found to provide the best 
match between the hybrids and the phenomenological 
model. 



B. Amplitude model 

In a similar manner to the phase, we approach the 
problem of fitting the amplitude of the GW by noting 
that the PN amplitude obtained from the SPA expression 
could be formally re-expanded as 



Ai^^if) = cn-y^ 1 



5 
fe=2 



n^/^ 



(5.10) 



where fl = nAIf. We introduce a higher-order term to 
model the pre-merger amplitude Apuif) 



Apm(/)=v4pn(/)+7i/' 



5/6 



(5.11) 



where the 71 coefficient is introduced to model the ampli- 
tude in the pre-merger regime and ApN is the amplitude 
constructed in Sec. Ill (see Fig. [2]). 

The ansatz for the amplitude during the ringdown is 

iRD(/) = ^i/:(/,/RD(a, M),,520(a)) /"'/', (5.12) 

where only the width and overall magnitude of the 
Lorentzian function £(/, /o, a) = a'^/ ((/ — /o)^ + o-^/4) 
are fitted to the hybrid data. The factor /^^^^ is in- 
troduced to correct the Lorentzian at high frequencies, 
since the hybrid data shows a faster fall-off, and Si ac- 
counts for the overall amplitude scale of the ringdown. In 
principle, the phenomenological parameter 62 should not 
be necessary because the width of the Lorentzian for the 
ringdown should be given by the quality factor Q which 
depends only on the spin of the final black hole. However, 
recall that here we estimate the final spin from the initial 
configuration using the fit given in [57]; 62 accounts for 
the errors in this fit. 

The phenomenological amplitude is constructed from 
these two pieces in a manner analogous to the phase 

4hcn(/) - ApuiDwJ^ + AnBif)w+, (5.13) 

with /o — 0.98/r,d and d = 0.015. Fig. |8] demonstrates 
how this phenomenological ansatz fits the hybrid am- 
plitude in a smooth manner through the late inspiral, 
merger and ringdown. 
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FIG. 9. Map of the phenomenological parameters to the physical parameters of the binary. 
TABLE II. CoefScionts to map the 9 free parameters of our phenomenological model to the physical parameters of the binary. 



Afe 


^(01) 


a(02) 


(^(11) 


^(10) 


^(20) 


Qfl 


-2.417 X 10^'^ 


-1.093 X 10^3 


-1.917 X 10^2 


7.267 X 10^2 


-2.504 X 10"' 


02 


5.962 X 10"^ 


-5.6 X 10"^ 


1.52 X 10~' 


-2.97 


1.312 X 10' 


03 


-3.283 X 10' 


8.859 


2.931 X 10' 


7.954 X 10' 


-4.349 X 10^ 


OtA 


1.619 X 10^ 


-4.702 X 10' 


-1.751 X 10^ 


-3.225 X 10^ 


1.587 X 10^ 


as 


-6.32 X 10^ 


2.463 X 10^ 


1.048 X 10^ 


3.355 X 10^ 


-5.115 X 10^ 


an 


-4.809 X 10' 


-3.643 X 10^ 


-5.215 X 10^ 


1.87 X 10^ 


7.354 X 10^ 


71 


4.149 


-4.07 


-8.752 X 10' 


-4.897 X 10' 


6.665 X 10^ 


5i 


-5.472 X 10"^ 


2.094 X 10"^ 


3.554 X 10-' 


1.151 X 10"' 


9.64 X 10"' 


h 


-1.235 


3.423 X 10"' 


6.062 


5.949 


-1.069 X 10' 



C. Mapping the phenomenological coefficients 

Our models for the amplitude and phase involve 9 phe- 
nomenological param eters {ai, ■ . . , ag, 71, (5i, (52} defined 
in Eqs. $^, Ku\ and (|5.12|). The coefficients /3i,2 



from (5.7) can be trivially derived from the set of a^- 
We now need to find the mapping A^ — > A^ from the 
physical to these phenomenological parameters. As men- 
tioned earlier, instead of xi,2, we co nsid er x, the weighted 
sum of the spins, defined in Eq. (5.1). Thus, our phe- 



nomenological waveforms are parametrized only by the 
symmetric mass ratio 77 and the spin parameter x, as 
well as by the total mass of the system M through a 
trivial rescaling. Fig. [9] shows the mapping of ak, "fk and 
(5fc to surfaces in the (77, x)-plane. 

The 9 phenomenological coefficients introduced in our 
model, denoted generically by A^ , are expressed in terms 
of the physical parameters of the binary as 



A. 



E 



cf^^^x^ 



(5.14) 



i+je{i,2} 



which yields 5 coefficients C'*-'' 
ters, as given in table [fT) 



for each of the 9 parame- 



We evaluate the goodness of fit between the phe- 
nomenological model and the hybrid waveforms in terms 
of the fitting factor, i.e. the ambiguity function A{X, A') 
defined in Eq. (4.10) and the overlap, i.e. O = ^(A, A). 



In evaluating the overlap, we maximize over the extrinsic 
parameters io:0o as indicated in Eq. (4.10), but for the 

we do not 



10 



results shown in the upper panel of Fig. 
perform the additional maximization over the model pa- 
rameters A'. Thus, the results shown there can be viewed 
as a lower bound on the effectualness. The maximization 
over the intrinsic parameters rj, x and M allows to study 
the faithfulness of the model. 

Figs. [To] and [Tl] illustrate the result using the design 
curve of the Advanced LIGO detector. Fig. [TO] shows 
the overlap and fitting factor between the hybrid wave- 
forms constructed in Sec. IV D and their corresponding 



phenomenological fit. The match approaches unity by 
construction at low masses and degrades with increasing 
total mass. Nevertheless, for none of the hybrid wave- 
forms employed in the construction of our model does the 
overlap fall below a value of ~ 0.97, thus reflecting the 
fact that the phenomenological model effectually repre- 
sents the target signals. A further maximization over the 
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FIG. 10. Overlaps and fitting factors between the hybrid 
waveform constructed according to the procedure described 
in Sec. |IVD| and the proposed phenomenological fit, using 
the design sensitivity curve of Advanced LIGO. The labels 
indicate the values of (77, x) for some configurations. In the 
upper panel We plot 0{X) = .4(A,A), i.e. we compute the 
ambiguity function (4.101 without maximizing over the pa- 
rameters of the model waveform; this is a lower bound on 
the effectualness. The bottom panel shows the maximized 
overlaps, i.e. .4(A, A'); the maximum bias of the optimized A' 
parameters is A77 = 5 x 10"^, Ax = 5 x 10"^, AM = 3 Mq. 



X' parameters, shown in the lower panel of Fig. 10 indi- 
cates a maximum bias on the intrinsic parameters of the 
binary of A77 = 5 x 10-^ Ax = 5 x IQ-"^, AM = 3 Mq. 

We have constructed a gravitational waveform model 
for binary black hole inspiral and coalescence starting 
with a particular set of simulations, and using a par- 
ticular ansatz for the waveform. Is this model robust, 
and is it consistent with waveforms from other numeri- 



cal simulations? In the upper panel of Fig. 11 and as 
a further test to assess the robustness of our model, we 
compute the maximized overlap between the phenomeno- 
logical waveforms and the NR data-sets #4-7 that were 
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FIG. 11. Upper panel: maximized overlaps ^(A, A') between 
the NR data-sets #4-7ab and the predicted phenomenological 
waveforms from our model for Advanced LIGO. The labels 
indicate the values of (77, x) for some configurations. Note that 
the short duration of the NR data prevents us from computing 
overlaps at lower masses. The maximization in this case has 
been done over 77 and x keeping M fixed and the maximum 
bias on the maximized parameters is Ar7 = 6 x 10"'', Ax = 
5 X 10"^ 



not used in the construction of the model. At low masses, 
there is no contribution of these short NR waveforms in 
the frequency band of interest for Advanced LIGO, and 
it turns out that the overlaps can be computed only for 
M> lOOAf©. 

In Fig. [11] we see that the maximization of the over- 
laps with respect to 77 and x shows values > 0.97 for 
all configurations; in this case the maximum bias in the 
parameters is A77 « 6 x 10~^,Ax « 5 x 10~^. This is 
roughly consistent with Fig. [TO] which shows the overlap 
and fitting factor of the model with the original set of 
hybrid waveforms. These results prove that our model 
is effectual and, thus, sufficient for detection. We shall 
study its effectualness and faithfulness in greater detail 
in a forthcoming paper. 



VI. SUMMARY AND FUTURE WORK 

The aim of this paper has been to construct an ana- 
lytical model for the inspiral and coalescence of binary 
black hole systems with aligned spins and comparable 
masses in circular orbits. Since this requires merging 
post-Newtonian and numerical relativity waveforms, one 
of the main themes has been to quantify the internal con- 
sistency of hybrid waveforms. This is important because 
even if one succeeds in finding a useful fit for a family of 
hybrid waveforms, one still needs to show that the hybrid 
one started with is a sufficiently good approximation to 
the true physical waveforms. We investigated the sys- 
tematics of constructing hybrid waveforms for accurate 
non-spinning waveforms based on the Llama code and 
we saw that neither the numerical errors nor the hybrid- 
construction errors are significant. This suggests that in 
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order to improve the accuracy of hybrid waveforms, we 
require either longer NR waveforms so that the matching 
with PN can be done earher in the inspiral phase, or im- 
proved PN models that are more accurate at frequencies 
closer to the binary merger. 

With the hybrid waveforms for non-precessing sys- 
tems in hand, we constructed an analytical model for 
the waveform which has an overlap and fitting factor 
of better than 97% for Advanced LIGO with the hy- 
brid waveforms for systems with a total mass ranging 
up to ~ 350Mq. Since these overlaps are comparable 
to those achieved with the alternative phenomenologi- 
cal waveform construction presented in |40) , we conclude 
that this process is robust, and, in particular, its accu- 
racy is not affected by the way in which the transitions 
between inspiral, merger and ringdown are modeled. Fur- 
thermore, though we have not discussed it in detail in this 
paper, it turns out that the model presented here agrees 
very well with the model of [JD]. This will be discussed 
in detail in a forthcoming paper |112j 

In the future we will study in greater detail the effectu- 
alness and faithfulness of this waveform model, thereby 
quantifying more precisely its performance for detection 
and parameter estimation. In this context it is impor- 
tant to extend this work to modes higher than the dom- 
inant £ ^ 2, m ^ ±2 spherical harmonics. It was shown 
recently |116| that the overlap with the real signal can 
possibly be affected by the inclusion of higher modes up 
to the order of ^ 1%, which is comparable or greater 
than the disagreement we find between hybrid and phe- 
nomenological model. 

We will further quantify the behavior of our templates 
in real non-Gaussian detector noise, and use them in 
real searches for gravitational wave signals. Eventu- 
ally, work is underway in extending the model to in- 



clude precessing spins. Our phenomenological model 
can be readily applied to existent GW detection efforts 
within the LIGO/Virgo Scientific Gollaborations. Ongo- 
ing searches are already making use of inspiral-merger- 
ringdown waveforms, such as the EOBNR family and the 
phenomenological family of [37ti40] in the form of soft- 
ware injections and as filter approximants. Our newly 
developed frequency-domain matching procedure should 
serve to cross-check the validity of these alternative ap- 
proaches and to complement them. 
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Appendix A: PN expansion coefficients 



For the convenience of the reader we explicitly give all the PN expansion coefficients used in Section [III] as functions 
of the symmetric mass ratio rj (3.2), the dimensionless spin magnitudes Xi = (•S'i • L)/mf, where L is the unit angular 



momentum vector, and x = Xi ttii/M -t- X2 milM. The energy (3.3) is given in terms of 
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■^E ~ 0.5772 is the Euler constant. Note that the next-to-leading order spin-orbit effects appearing at relative 2.5PN 
order (/s) have recently been corrected [117] and we take these corrections into account. 
The Taylor T4 approximant can be written as a series (3.6 1 with the following coefhcients 
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The spin-dependent terms that appear at 3 and 3.5 PN order (i.e., in ag and ay) are not complete since the corre- 
sponding terms are not known in energy and flux. However, in this re-expansion they do appear as contributions from 
lower order spin effects and we keep them. 



The TaylorF2 description of the Fourier phase (3.13) is expressed in terms of 
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The comment we just made about the spin contributions at 3 and 3.5 PN order holds for the a-coefficients of the 
TaylorF2 phase as weh. Also, note that the contributions in as that are not proportional to ln(7r/) could be absorbed 
in a re-definition of the undetermined additional phase (f>o that appears in Eq. (3.13). (A similar discussion can be 
found in 75 .) However, since we chose to set (I)q — when combining this phase description with other analytical 
formulas in our phenomenological model ( |5.9[ ), it is important to keep all terms in a^. 
The time-domain amplitude coefficients collected from [261 ESI EZ] read 
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